6 - Cell-Cell Communication analysis

Author

CDN team

Published

March 8, 2024

Introduction

In this notebooks we are going to carry out a cell-cell communication analysis using CellChat. You can see the GitHub repository here, the CellChat vignette here and the differential communication analysis vignette here. The goal of the notebook is to identify which communication pathways are altered between flu and covid infection!

Some associated literature which is a must read are:

Key Takeaways

  • The Ligand-Receptor (L-R) database used gather the collective prior knowledge and has a great impact on the results obtained.

  • Different CCC tools have varying assumptions, therefore, the tool of choise will also have a major impact on the results.

    • CellChat and CellPhoneDB make a point of modelling CCC events taking into account heteromeric complexes. This ensures all the subunits of a protein complex are expressed to consider a cell-cell interaction feasible. This assumption reduces false positive predictions.

    • CellChat, additionally, accounts for interaction mediator proteins such as agonists.

  • Broadly, CCC tools are generally able to capture relevant biological signals. However, predicted interactions tend to have false positives, if available leveraging information from additional modalities and analyses could help to refine the predictions.

  • CCC inference from scRNAseq data makes the assumption that gene expression is a proxy for protein levels. Moreover, they don’t (and can’t) account for other intermediate steps between translation and protein function such as post-translational modifications, secretion, diffusion…

Library

options(future.globals.maxSize= 891289600)
### Make sure all the packages are installed
if (!requireNamespace("Seurat", quietly = TRUE))
    install.packages("Seurat")

if (!requireNamespace("dplyr", quietly = TRUE))
    install.packages("dplyr")

if (!requireNamespace("devtools", quietly = TRUE))
    install.packages("devtools")
devtools::install_github("jinworks/CellChat")

if (!requireNamespace("ComplexHeatmap", quietly = TRUE))
    BiocManager::install("ComplexHeatmap")

### Load all the necessary libraries
library(Seurat)
library(dplyr)
library(CellChat)
library(ComplexHeatmap)

Load data

We’re going to be working with a dataset from the paper - Immunophenotyping of COVID-19 and influenza highlights the role of type I interferons in development of severe COVID-19 Download data from cellxgene portal.

# Download the data in data/ directory
# download.file(
#     url = "https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds",
#     destfile = "../data/workshop-data.rds",
#     method = "wget",
#     extra = "-r -p --random-wait")
# We can also use the CLI with the wget command below
# wget https://datasets.cellxgene.cziscience.com/d8e35450-de43-451a-9979-276eac688bce.rds

se <- readRDS("../data/workshop-data.rds")

Analysis

Convert ENSEMBL IDs to Gene Symbols

Right away we can see how ensembl ids are used in the rownames. Let’s transform them into their matched symbols to make them human-readable:

head(rownames(se))
[1] "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457"
[5] "ENSG00000000460" "ENSG00000000938"

Convert to gene symbols

gene_df <- readr::read_csv(file = "../data/cov_flu_gene_names_table.csv")

symbol_id <- data.frame(index = rownames(se)) %>%
    left_join(gene_df) %>%
    pull(feature_name)

# re-create seurat object
mtx <- se@assays$RNA@data
rownames(mtx) <- symbol_id
se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
# symbol_id <- mapIds(
#     org.Hs.eg.db,
#     keys = rownames(se), 
#     column = "SYMBOL",
#     keytype = "ENSEMBL",
#     multiVals = "first")
# 
# # df <- data.frame(symbol = symbol_id, ensembl = names(symbol_id))
# all(rownames(se) == names(symbol_id))
# 
# # re-create seurat object
# mtx <- se@assays$RNA@data
# rownames(mtx) <- symbol_id
# 
# # Remove NAs
# sum(is.na(rownames(mtx)))
# dim(mtx)
# mtx <- mtx[!is.na(rownames(mtx)), ]
# dim(mtx)
# 
# rownames(mtx) <- make.unique(rownames(mtx))
# se <- CreateSeuratObject(counts = mtx, meta.data = se@meta.data)
# 
# sum(is.na(rownames(se)))

Split dataset

# First we will remove the Uncategorized1/2 populations as we have seen that they are doublets
se <- se[, ! se$Celltype %in% c("Uncategorized1", "Uncategorized2")]

# Create samples so the CellChat knows which cells come from each sample
se$samples <- se$`Sample ID`
# Convert to character to drop inexisting factors
se$Celltype <- factor(as.character(se$Celltype))
# Subset to cells coming form COVID-19 patients
covid <- se[, se$disease == "COVID-19"]

# Subset to cells coming form flu patients
flu <- se[, se$disease == "influenza"]

Let’s visualize the cell type populations between these 2 datasets:

bind_rows(
    data.frame(prop.table(table(covid$Celltype)), disease = "covid"),
    data.frame(prop.table(table(flu$Celltype)), disease = "flu")) %>%
    ggplot(aes(x = Var1, y = disease, fill = Freq, label = round(Freq, 2))) +
    geom_tile(color = "lightgrey") +
    geom_text(color = "lightgrey") +
    scale_fill_viridis_c(option = "magma") +
    labs(title = "Cell type proportions", x = "", y = "") +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

Carry out pre-processing steps on covid data and create the cellchat object

# Normalize data
covid <- NormalizeData(covid, verbose = FALSE)

cc_covid <- createCellChat(
    object = covid,
    group.by = "Celltype",
    assay = "RNA",
    do.sparse = TRUE)
[1] "Create a CellChat object from a Seurat object"
The `meta.data` slot in the Seurat object is used as cell meta information 
Set cell identities for the new CellChat object 
The cell groups used for CellChat analysis are  B cell, IgG-, B cell, IgG+, CD4, EM-like, CD4, non-EM-like, CD8, EM-like, CD8, non-EM-like, classical Monocyte, DC, intermediate Monocyte, NK cell, nonclassical Monocyte, Platelet, RBC 

Same for the flu

# Normalize data
flu <- NormalizeData(flu, verbose = FALSE)

cc_flu <- createCellChat(
    object = flu,
    group.by = "Celltype",
    assay = "RNA",
    do.sparse = TRUE)
[1] "Create a CellChat object from a Seurat object"
The `meta.data` slot in the Seurat object is used as cell meta information 
Set cell identities for the new CellChat object 
The cell groups used for CellChat analysis are  B cell, IgG-, B cell, IgG+, CD4, EM-like, CD4, non-EM-like, CD8, EM-like, CD8, non-EM-like, classical Monocyte, DC, intermediate Monocyte, NK cell, nonclassical Monocyte, Platelet, RBC 

CellChat CCC analysis

Cellchat intuition

How is the communication probability computed? See the methods section in the CellChat paper.

The communication probability Pi,j from cell groups i to j for a particular ligand-receptor pair k was modeled by

\[ P_{i,j}^k = \frac{{L_iR_j}}{{K_h + L_iR_j}} \times \left( {1 + \frac{{AG_i}}{{K_h + AG_i}}} \right) \cdot \left( {1 + \frac{{AG_j}}{{K_h + AG_j}}} \right) \\ \times \frac{{K_h}}{{K_h + AN_i}} \cdot \frac{{K_h}}{{K_h + AN_j}} \times \frac{{n_in_j}}{{n^2}} \]

\[ {L_i = \root {{m1}} \of {{L_{i,1} \cdots L_{i,m1}}},\,R_j = \root {{m2}} \of {{R_{j,1} \cdots R_{j,m2}}} \cdot \frac{{1 + RA_j}}{{1 + RI_j}}} \]

  • Li and Rj represent the expression level of ligand L and receptor R in cell group i and cell group j, respectively

  • Kh whose default value was set to be 0.5

  • AG is the average expression of the L-R agonists

CellChat Database

Define the use of the human L-R database. The database used is a key decision in the analysis workflow since it will define the universe of Ligand-Receptor interactions we are going to be looking for. We want it to be as comprehensive as possible while simultaneously only including curated communication partners.

We personally like the CellChatDB database because:

  • Groups LR interactions into functional pathways.

  • Takes into account multimeric complexes as well as other cofactors, such as agonists and antagonists present in the dataset.

  • LR interactions come from KEGG Pathway database and manually curated from primary references.

We can take a glimpse at all of this information by looking at the database loaded into R!

CellChatDB <- CellChatDB.human
showDatabaseCategory(CellChatDB)

# glimpse(CellChatDB)
CellChatDB$interaction %>%
    dplyr::select(interaction_name, pathway_name, ligand, receptor, agonist, evidence, annotation) %>% 
    head() %>%
    as_tibble()
# A tibble: 6 × 7
  interaction_name    pathway_name ligand receptor   agonist evidence annotation
  <chr>               <chr>        <chr>  <chr>      <chr>   <chr>    <chr>     
1 TGFB1_TGFBR1_TGFBR2 TGFb         TGFB1  TGFbR1_R2  TGFb a… KEGG: h… Secreted …
2 TGFB2_TGFBR1_TGFBR2 TGFb         TGFB2  TGFbR1_R2  TGFb a… KEGG: h… Secreted …
3 TGFB3_TGFBR1_TGFBR2 TGFb         TGFB3  TGFbR1_R2  TGFb a… KEGG: h… Secreted …
4 TGFB1_ACVR1B_TGFBR2 TGFb         TGFB1  ACVR1B_TG… TGFb a… PMID: 2… Secreted …
5 TGFB1_ACVR1C_TGFBR2 TGFb         TGFB1  ACVR1C_TG… TGFb a… PMID: 2… Secreted …
6 TGFB2_ACVR1B_TGFBR2 TGFb         TGFB2  ACVR1B_TG… TGFb a… PMID: 2… Secreted …
CellChatDB$interaction %>%
    dplyr::select(interaction_name, pathway_name, ligand, receptor, agonist, evidence, annotation) %>% 
    tail() %>%
    as_tibble()
# A tibble: 6 × 7
  interaction_name      pathway_name ligand receptor agonist evidence annotation
  <chr>                 <chr>        <chr>  <chr>    <chr>   <chr>    <chr>     
1 Testosterone-Testost… Testosterone Testo… AR       ""      HMRbase… Non-prote…
2 Testosterone-Testost… Testosterone Testo… AR       ""      HMRbase… Non-prote…
3 Thyrostimulin-GPHB5_… TSH          GPHB5… TSHR     ""      PMID: 1… Non-prote…
4 Thyrotropin-CGA_TSHB… TSH          CGA_T… TSHR     ""      PMID: 1… Non-prote…
5 Triiodothyronine-T3-… Triiodothyr… T3-DI… THRA     ""      HMRbase… Non-prote…
6 Triiodothyronine-T3-… Triiodothyr… T3-DI… THRB     ""      HMRbase… Non-prote…

Preprocess expression data

To infer the cell state-specific communications, CellChat identifies over-expressed ligands or receptors in one cell group and then identifies over-expressed ligand-receptor interactions if either ligand or receptor are over-expressed.

# subset the expression data of signaling genes for saving computation cost
cc_covid@DB <- CellChatDB
cc_covid <- subsetData(cc_covid) # This step is necessary even if using the whole database
future::plan("multisession", workers = 4) # do parallel
cc_covid <- identifyOverExpressedGenes(cc_covid)
cc_covid <- identifyOverExpressedInteractions(cc_covid)
The number of highly variable ligand-receptor pairs used for signaling inference is 2205 
# The number of highly variable ligand-receptor pairs used for signaling inference is 2212 

Same for the flu

# subset the expression data of signaling genes for saving computation cost
cc_flu@DB <- CellChatDB
cc_flu <- subsetData(cc_flu) # This step is necessary even if using the whole database
future::plan("multisession", workers = 4) # do parallel
cc_flu <- identifyOverExpressedGenes(cc_flu)
cc_flu <- identifyOverExpressedInteractions(cc_flu)
The number of highly variable ligand-receptor pairs used for signaling inference is 1720 
# The number of highly variable ligand-receptor pairs used for signaling inference is 1728 

Inference of cell-cell communication network

From the CellChat vignette: CellChat infers the biologically significant cell-cell communication by assigning each interaction with a probability value and peforming a permutation test. CellChat models the probability of cell-cell communication by integrating gene expression with prior known knowledge of the interactions between signaling ligands, receptors and their cofactors using the law of mass action.

The number of inferred ligand-receptor pairs clearly depends on the method for calculating the average gene expression per cell group. trimean approximates 25% truncated mean, implying that the average gene expression is zero if the percent of expressed cells in one group is less than 25%.

By setting this stringent parameters we ensure we are capturing robust signal and removing spurious cell-cell communication.

cc_covid <- computeCommunProb(cc_covid, type = "triMean")
triMean is used for calculating the average gene expression per cell group. 
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2024-03-08 10:32:05.470454]"
[1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2024-03-08 10:34:34.199726]"
cc_flu <- computeCommunProb(cc_flu, type = "triMean")
triMean is used for calculating the average gene expression per cell group. 
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2024-03-08 10:34:36.527255]"
[1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2024-03-08 10:36:04.791565]"
# If we wanted to be more lax we could use truncatedMean
# By setting trim - 0.1 we are only removing the top and bottom 10% expressing cells for that gene
# cc_covid <- computeCommunProb(cc_covid, type = "truncatedMean", trim = 0.1)
# cc_flu <- computeCommunProb(cc_flu, type = "truncatedMean", trim = 0.1)

Summarize CCC to a signaling pathway level

We can summarize the individual ligand-receptor interactions into pathways for a broader and more functional look at the data.

cc_covid <- computeCommunProbPathway(cc_covid)
cc_flu <- computeCommunProbPathway(cc_flu)

## Look at all the pathways availavle
cc_covid@netP$pathways
 [1] "MHC-I"         "MHC-II"        "MIF"           "CD99"         
 [5] "GALECTIN"      "CypA"          "CLEC"          "ICAM"         
 [9] "ANNEXIN"       "ADGRE"         "CD45"          "CCL"          
[13] "APP"           "Prostaglandin" "BAFF"          "SELPLG"       
[17] "LCK"           "PECAM1"        "CD23"          "RESISTIN"     
[21] "TGFb"          "Cholesterol"   "VISFATIN"      "IL16"         
[25] "SEMA4"         "PARs"          "GRN"           "TNF"          
[29] "CXCL"          "GP1BA"         "COLLAGEN"      "CD40"         
[33] "CysLTs"        "CD160"         "CD86"          "LAIR1"        
[37] "CD6"           "PECAM2"        "IFN-II"        "ESAM"         
[41] "ICOS"          "BAG"           "FLT3"          "IL1"          
[45] "TWEAK"         "NECTIN"        "MPZ"          

Calculate the aggregated cell-cell communication network by counting the number of links or summarizing the communication probability.

cc_covid <- aggregateNet(cc_covid)
cc_flu <- aggregateNet(cc_flu)

Lastly we can compute the network centrality scores of the cell types which we can use to identify key senders, receivers, mediators…

# the slot 'netP' means the inferred intercellular communication network of signaling pathways
cc_covid <- netAnalysis_computeCentrality(cc_covid, slot.name = "netP") 
cc_flu <- netAnalysis_computeCentrality(cc_flu, slot.name = "netP") 

Visualize inferred IFN-II signaling network

With Heatmap

object.list <- list(covid = cc_covid, flu = cc_flu)

par(mfrow = c(1, 2))
plt_ls <- lapply(seq_len(length(object.list)), function(i) {
  netVisual_heatmap(
    object.list[[i]],
    signaling = "IFN-II",
    color.heatmap = "Reds",
    title.name = glue::glue("IFN-II signaling network - {names(object.list)[i]}"))
})

draw(plt_ls[[1]] + plt_ls[[2]], ht_gap = unit(.5, "cm"))

Merge Covid and Flue objects

Now we can merge the CellChat objects and work with them together. The rest of the analysis follows the comparison vignette.

cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)
The cell barcodes in merged 'meta' is  AAACCCAAGAATGTTG-12 AAACCCAAGGCCTGCT-12 AAACCCAAGGGCAATC-1 AAACCCAAGTAAGGGA-18 AAACCCACAAGAGCTG-17 AAACCCACACAGTCGC-20 

Save objects in case we need them for later use

saveRDS(object.list, file = "../data/cellchat_ls.rds")
# object.list <- readRDS(file = "../data/cellchat_ls.rds")
# rm(object.list); gc()

Differential interaction analysis

Let’s start by comparing the the total number of interactions and the sum of the weights

gg1 <- compareInteractions(
    cellchat,
    show.legend = FALSE,
    group = c(1,2))
gg2 <- compareInteractions(
    cellchat,
    show.legend = FALSE,
    group = c(1,2), 
    measure = "weight")
gg1 + gg2

With this plot we can visualize how the covid dataset has a slightly larger number of interactions but the interaction strength within the flu dataset are stronger/more likely.

We can also take a look at how these interactions are different between covid and flu. In the colorbar, red (or blue) represents increased (or decreased) signaling in the second dataset compared to the first one. In our case red is stronger in Flu and blue in covid

gg1 <- netVisual_heatmap(cellchat)
gg2 <- netVisual_heatmap(cellchat, measure = "weight")

gg1 + gg2

Next we can see how much the incoming and outgoing interaction signals are changing between conditions. This requires some data prep so we are showing the plot below:

slot.name <- "netP"
x.measure <- "outdeg"
y.measure <- "indeg"
signaling <- NULL

df_ls <- lapply(names(object.list), function(nm) {
    object <- object.list[[nm]]
    centr <- slot(object, slot.name)$centr
    outgoing <- matrix(0, nrow = nlevels(object@idents), ncol = length(centr))
    incoming <- matrix(0, nrow = nlevels(object@idents), ncol = length(centr))
    dimnames(outgoing) <- list(levels(object@idents), names(centr))
    dimnames(incoming) <- dimnames(outgoing)
    for (i in 1:length(centr)) {
        outgoing[, i] <- centr[[i]][[x.measure]]
        incoming[, i] <- centr[[i]][[y.measure]]
    }
    outgoing.cells <- rowSums(outgoing)
    incoming.cells <- rowSums(incoming)
    num.link <- aggregateNet(object, signaling = signaling, return.object = FALSE, 
                             remove.isolate = FALSE)$count
    num.link <- rowSums(num.link) + colSums(num.link) - diag(num.link)
    df <- data.frame(
        outgoing.cells,
        incoming.cells, 
        names(incoming.cells),
        num.link,
        nm)
    colnames(df) <- c(paste0("outgoing_", nm), paste0("incoming_", nm), "labels", paste0("count_", nm))
    df$labels <- factor(df$labels, levels = names(incoming.cells))
    
    df
})

Visualize how the interactions changed

require(gridExtra)
library(colorBlindness)
library(RColorBrewer)
# Define the number of colors you want
nb.cols <- length(unique(se$Celltype))
col_pal <- colorRampPalette(brewer.pal(8, "Set2"))(nb.cols)
names(col_pal) <- sort(unique(se$Celltype))
# Create a ggplot with 18 colors 
df_ls[[1]] %>%
    left_join(df_ls[[2]], by = "labels") %>%
    
    ggplot() +
    geom_point(aes(x = outgoing_flu, y = incoming_flu, color = labels,
                   fill = "Flu", size = count_flu), alpha = 0.5) +
    geom_point(aes(x = outgoing_covid, y = incoming_covid, color = labels,
                   fill = "Covid", size = count_covid), size = 3, pch = 21) +
    geom_segment(aes(
        x = outgoing_covid,
        y = incoming_covid,
        xend = outgoing_flu,
        yend = incoming_flu,
        color = labels),
        arrow = arrow(angle = 40, length = grid::unit(.35, "cm")),
        show.legend = FALSE) +
    theme_classic() +
    guides(
        size = FALSE,
        shape = FALSE,
        colour = guide_legend(override.aes = list(size = 5))) +
    scale_color_manual(values = col_pal) +
    labs(
        title = "Cell-Cell interaction strength shifts from Covid to Flu",
       # x = "Outgoing strength",
      #  y = "Incoming strength",
      y = NULL, x = NULL,
        color = "Cell Type",
        fill = "Disease")  +
        scale_fill_manual(values = c("white", "black")) +
    theme(panel.grid.minor = element_blank(),
          panel.grid.major = element_line()) +
    scale_y_continuous(labels = function(x) {
        case_when(x == 25 ~ "25 Incoming\nStrength",
                  TRUE ~ as.character(round(x)))
    })  +
    scale_x_continuous(labels = function(x) {
        case_when(x == 10 ~ "10\nOutgoing\nStrength",
                  TRUE ~ as.character(round(x)))
    }) +
    coord_fixed()

In this case we can see how CD8 T cells greatly increase their incoming signal as well as intermediate monocytes. NKs and CD4 EM-like, in turn, greatly increase their outgoing signal while not increasing their incoming signal.

CellChat allows us to explore signalling changes at the pathway level within one specific population between both conditions. Positive values highlight increases in the second group (Flu) while negative values are increases in the first group (covid).

netAnalysis_signalingChanges_scatter(cellchat, idents.use = "CD4, EM-like")

Differential pathway signalling

Look at broad changes in pathway signalling across the entire dataset

gg1 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = F, do.stat = TRUE)

gg1 + gg2

Diving into a pathway

We also need to take a look at which are the LR interactions driving a pathway. We can easily look at this using the function XX

(netAnalysis_contribution(object.list[[1]], signaling = "MHC-II") +
    labs(title = "Contributions of each L-R pair in Covid")) / 
(netAnalysis_contribution(object.list[[2]], signaling = "MHC-II") +
    labs(title = "Contributions of each L-R pair in Flu"))

Comparing all the pathways at once

Lastly, we can also take a look at all the pathway changes at once by cell type between both conditions.

library(ComplexHeatmap)

# combining all the identified signaling pathways from different datasets 
pathway.union <- union(object.list[["covid"]]@netP$pathways, object.list[["flu"]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(
    object.list[["covid"]],
    pattern = "outgoing",
    signaling = pathway.union,
    title = "covid",
    width = 5,
    height = 12)

ht2 <- netAnalysis_signalingRole_heatmap(
    object.list[["flu"]],
    pattern = "outgoing",
    signaling = pathway.union,
    title = "flu",
    width = 5,
    height = 12)

draw(ht1 + ht2, ht_gap = unit(.5, "cm"))

# combining all the identified signaling pathways from different datasets 
pathway.union <- union(object.list[["covid"]]@netP$pathways, object.list[["flu"]]@netP$pathways)
ht1 <- netAnalysis_signalingRole_heatmap(
    object.list[["covid"]],
    pattern = "incoming",
    signaling = pathway.union,
    title = "covid",
    width = 5,
    height = 12)

ht2 <- netAnalysis_signalingRole_heatmap(
    object.list[["flu"]],
    pattern = "incoming",
    signaling = pathway.union,
    title = "flu",
    width = 5,
    height = 12)

draw(ht1 + ht2, ht_gap = unit(.5, "cm"))

Visualizing L-R interactions between cell types

levels(cellchat@idents$joint)
 [1] "B cell, IgG-"          "B cell, IgG+"          "CD4, EM-like"         
 [4] "CD4, non-EM-like"      "CD8, EM-like"          "CD8, non-EM-like"     
 [7] "classical Monocyte"    "DC"                    "intermediate Monocyte"
[10] "NK cell"               "nonclassical Monocyte" "Platelet"             
[13] "RBC"                  
netVisual_bubble(
    cellchat,
    sources.use = "DC",
    targets.use = c("NK cell", "CD4, EM-like"),
    comparison = c(1, 2),
    angle.x = 45)

We can see here how the MHC-II signalling is present between DCs & CD4 and completely absent with NKs.

Session Info

sessionInfo()
R version 4.3.1 (2023-06-16)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/New_York
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] RColorBrewer_1.1-3    colorBlindness_0.1.9  gridExtra_2.3        
 [4] ComplexHeatmap_2.16.0 CellChat_2.1.2        Biobase_2.60.0       
 [7] BiocGenerics_0.46.0   ggplot2_3.5.0         igraph_2.0.2         
[10] dplyr_1.1.4           Seurat_5.0.2          SeuratObject_5.0.1   
[13] sp_2.1-3             

loaded via a namespace (and not attached):
  [1] RcppAnnoy_0.0.22       splines_4.3.1          later_1.3.2           
  [4] tibble_3.2.1           polyclip_1.10-6        ggnetwork_0.5.13      
  [7] fastDummies_1.7.3      lifecycle_1.0.4        rstatix_0.7.2         
 [10] doParallel_1.0.17      vroom_1.6.3            globals_0.16.2        
 [13] processx_3.8.2         lattice_0.21-8         MASS_7.3-60           
 [16] backports_1.4.1        magrittr_2.0.3         sass_0.4.8            
 [19] plotly_4.10.4          rmarkdown_2.25         jquerylib_0.1.4       
 [22] yaml_2.3.8             remotes_2.4.2.1        httpuv_1.6.14         
 [25] NMF_0.26               sctransform_0.4.1      spam_2.10-0           
 [28] sessioninfo_1.2.2      pkgbuild_1.4.2         spatstat.sparse_3.0-3 
 [31] reticulate_1.35.0.9000 cowplot_1.1.3          pbapply_1.7-2         
 [34] abind_1.4-5            pkgload_1.3.2.1        Rtsne_0.17            
 [37] presto_1.0.0           purrr_1.0.2            circlize_0.4.15       
 [40] IRanges_2.34.1         S4Vectors_0.38.1       ggrepel_0.9.5         
 [43] irlba_2.3.5.1          listenv_0.9.1          spatstat.utils_3.0-4  
 [46] goftest_1.2-3          RSpectra_0.16-1        spatstat.random_3.2-2 
 [49] fitdistrplus_1.1-11    parallelly_1.37.0      svglite_2.1.3         
 [52] leiden_0.4.3.1         codetools_0.2-19       tidyselect_1.2.0      
 [55] shape_1.4.6            farver_2.1.1           matrixStats_1.2.0     
 [58] stats4_4.3.1           spatstat.explore_3.2-6 jsonlite_1.8.8        
 [61] BiocNeighbors_1.18.0   GetoptLong_1.0.5       ellipsis_0.3.2        
 [64] progressr_0.14.0       ggalluvial_0.12.5      ggridges_0.5.6        
 [67] survival_3.5-7         iterators_1.0.14       systemfonts_1.0.4     
 [70] foreach_1.5.2          tools_4.3.1            sna_2.7-2             
 [73] ica_1.0-3              Rcpp_1.0.12            glue_1.7.0            
 [76] xfun_0.42              usethis_2.2.2          withr_3.0.0           
 [79] BiocManager_1.30.22    fastmap_1.1.1          fansi_1.0.6           
 [82] callr_3.7.3            digest_0.6.34          gridGraphics_0.5-1    
 [85] R6_2.5.1               mime_0.12              colorspace_2.1-0      
 [88] Cairo_1.6-1            scattermore_1.2        tensor_1.5            
 [91] spatstat.data_3.0-4    utf8_1.2.4             tidyr_1.3.1           
 [94] generics_0.1.3         data.table_1.15.0      FNN_1.1.4             
 [97] prettyunits_1.1.1      httr_1.4.7             htmlwidgets_1.6.4     
[100] uwot_0.1.16            pkgconfig_2.0.3        gtable_0.3.4          
[103] registry_0.5-1         lmtest_0.9-40          htmltools_0.5.7       
[106] carData_3.0-5          profvis_0.3.8          dotCall64_1.1-1       
[109] clue_0.3-64            tidyverse_2.0.0        scales_1.3.0          
[112] png_0.1-8              knitr_1.45             rstudioapi_0.15.0     
[115] tzdb_0.4.0             reshape2_1.4.4         rjson_0.2.21          
[118] coda_0.19-4.1          statnet.common_4.9.0   nlme_3.1-163          
[121] curl_5.2.0             cachem_1.0.8           zoo_1.8-12            
[124] GlobalOptions_0.1.2    stringr_1.5.1          KernSmooth_2.23-22    
[127] parallel_4.3.1         miniUI_0.1.1.1         pillar_1.9.0          
[130] vctrs_0.6.5            RANN_2.6.1             ggpubr_0.6.0          
[133] urlchecker_1.0.1       promises_1.2.1         car_3.1-2             
[136] xtable_1.8-4           cluster_2.1.4          evaluate_0.23         
[139] magick_2.7.5           readr_2.1.4            cli_3.6.2             
[142] compiler_4.3.1         rlang_1.1.3            crayon_1.5.2          
[145] rngtools_1.5.2         ggsignif_0.6.4         future.apply_1.11.1   
[148] labeling_0.4.3         ps_1.7.5               plyr_1.8.9            
[151] fs_1.6.3               stringi_1.8.3          network_1.18.2        
[154] BiocParallel_1.34.2    viridisLite_0.4.2      deldir_2.0-2          
[157] gridBase_0.4-7         munsell_0.5.0          lazyeval_0.2.2        
[160] devtools_2.4.5         spatstat.geom_3.2-8    Matrix_1.6-5          
[163] RcppHNSW_0.6.0         hms_1.1.3              patchwork_1.2.0       
[166] bit64_4.0.5            future_1.33.1          shiny_1.8.0           
[169] ROCR_1.0-11            broom_1.0.5            memoise_2.0.1         
[172] bslib_0.6.1            bit_4.0.5